Ce document présente un usage possible des packages R qui interfacent OpenStreetMap (OSM) : géocoder une adresse, import de tuiles et de couches géographiques, import de points d’intérêt, routage et calcul d’itinéraires.
Il propose aussi des méthodes spatiales (carroyages, KDE, segments de lignes) et des visualisations associées pour restituer graphiquement la structure spatiale de semis de points ou de lignes.
En utilisant exclusivement des informations issues d’OSM, l’exercice consiste à identifier les secteurs accessibles à pieds dans le voisinage géographique d’un lieu défini par une adresse.
Reproduire l’analyse
Pour reproduire l’analyse, télécharger l’archive ZIP associée au dépôt https://github.com/riatecom/osm-elementr-2025-application, ouvrir le fichier proj.Rproj puis le fichier script.R qui reproduit en tout point le contenu du support de formation.
Cette analyse est reproductible en n’importe quel lieu, avec une adresse (adr) et un nom de ville (lib). Néanmoins les temps de calcul utilisant le serveur de démonstration d’OSRM ont été paramétrés pour ne pas le surcharger. L’exécution de ces blocs de code prennent du temps.
En conséquence, 6 jeux de données ont été préparés en amont à partir de l’adresse des gares SNCF de Lille, Toulouse, Grenoble, Nantes, Clermont-Ferrand et Compiègne, ainsi que de la station de métro Front Populaire. Les blocs de code relevant des calculs d’itinéraires / temps de trajets sont placés en commentaires dans le fichier script.R, mais sont fonctionnels.
lib adr
1 Grenoble 1 Place de la Gare, 38000 Grenoble, France
2 Clermont 40 Avenue de l'Union Soviétique, 63000 Clermont-Ferrand, France
3 Toulouse Gare SNCF Matabiau, 31500 Toulouse, France
4 Compiegne Gare de Compiègne, 60200 Compiègne, France
5 Condorcet Front Populaire, 93210 Saint-Denis, France
6 Nantes 27 Boulevard de Stalingrad, 44000 Nantes, France
7 Lille Lille Flandres, 59000 Lille, France
Nous présentons ici les résultats depuis la gare SNCF de Compiègne.
# 1. Packages qui interfacent des données OSM# install.packages('maposm', repos = 'https://riatelab.r-universe.dev')library(tidygeocoder) # géocodagelibrary(maptiles) # import de tuiles OSM (raster)library(osmdata) # import de données OSM (vecteur)library(maposm) # import de données OSM (couches géo)library(osrm) # calcul d'itinéraires# 2. Autres packages utilitaireslibrary(sf) # manipulation de données vectorielleslibrary(terra) # manipulation de données rasterlibrary(mapsf) # cartographie thématiquelibrary(maplegend) # légendeslibrary(spatstat) # analyse de semis de pointslibrary(stplanr) # segmentiser des ligneslibrary(mapview) # cartographie interactivelibrary(lwgeom)
1. Géocodage
Le package tidygeocoder(Cambon et al. 2025) permet d’utiliser des services de géocodage en ligne. Par défaut, le package utilise Nominatim qui se base sur les données OSM.
loc <-"Mairie de Compiègne"pt <-data.frame(address ="28 Place de l'Hôtel de Ville, 60200 Compiègne")pt <-geocode(.tbl = pt, address ="address", quiet =TRUE) |>st_as_sf(coords =c("long", "lat"), crs =4326)pt
Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 2.826361 ymin: 49.41777 xmax: 2.826361 ymax: 49.41777
Geodetic CRS: WGS 84
# A tibble: 1 × 2
address geometry
* <chr> <POINT [°]>
1 28 Place de l'Hôtel de Ville, 60200 Compiègne (2.826361 49.41777)
2. Import de tuiles OSM
Le package maptiles(Giraud 2025) permet de télécharger des tuiles chez différents fournisseurs, dont OpenStreetMap, au format raster.20 niveaux de zoom sont disponibles, selon la taille de l’objet géographique à cartographier. Un zoom de niveau 10 correspond à une échelle 1/500.000.
Plusieurs modèles de tuiles sont disponibles ici, certains nécessitent une inscription.
# Autre exemple de fournisseurtiles <-get_tiles(x = emprise,project =FALSE,crop =TRUE,provider ="OpenStreetMap.HOT",zoom =14)mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt = loc)mf_title("Tuiles OSM - OpenStreetmap.HOT")
3. Import de couches géographiques OSM
Pour la cartographie, le package maposm (pas sur le CRAN) se base sur les fonctions du package osmdata(Padgham et al. 2023) (partie suivante) pour télécharger les données OSM, au format vectoriel.
La fonction principale, om_get(), permet de récupérer certaines couches géographiques utiles (routes, parcs, bâti, etc) autour d’un point, dans un rayon donné. Le résultat prend la forme d’une liste d’objets sf.
res <-om_get(x =c(st_coordinates(pt)[1],st_coordinates(pt)[2]),r =2000)
mf_map(res$zone, col ="#f2efe9", border =NA, add =FALSE)mf_map(res$urban, col ="#e0dfdf", border ="#e0dfdf", lwd = .5, add =TRUE)mf_map(res$green, col ="#c8facc", border ="#c8facc", lwd = .5, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$building, col ="#d9d0c9", border ="#c6bab1", lwd = .5, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Autour de la mairie")mf_arrow()mf_scale(size =500, scale_units ="m")mf_credits("OpenStreetMap contributors, 2025", pos ="bottomleft")
4. Export de tags avec osmdata
Le package osmdata permet d’extraire des données OSM à partir de l’API Overpass. Il nécessite dans un premier temps d’établir une bounding box au sein de laquelle seront extraites les données, puis de l’exécuter (fonction add_osm_feature()). Le résultat est une liste, qu’il faut ensuite convertir en objet sf.
Les objets OSM prennent la forme de nodes (points), ways (lignes ou polygones), ou relations (relation entre plusieurs objets). Ces éléments sont décrits par un ensemble de tags, recensées sous forme de clé-valeurs (key-values). Ici, nous recherchons dans la clé amenity les points (nodes) correspondant aux restaurants, bars et cafés. En dehors des clés-valeurs principales, d’autres attributs peuvent être présents, comme en l’occurence le type de cuisine, ou bien les horaires.
La fonction om_map() affiche les couches géographiques récupérées.
# Affichageom_map(res, title ="Autour de la mairie")mf_map(poi_reproj, type ="typo", var ="amenity", cex = .7, pch =15, leg_pos ="topright", leg_val_cex = .9, leg_title_cex = .9,add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, lwd =2, col ="darkblue",add =TRUE)
Il est possible d’utiliser le package mapview(Appelhans et al. 2023) pour afficher les points ainsi que leurs attributs, de manière interactive.
mapview(res$zone, alpha.regions =0) +mapview(poi)
osmextract
Le package osmextract(Gilardi et Lovelace 2024) permet d’extraire des données directement depuis une base de données OSM, sans avoir à passer par une API. Le plus petit niveau géographique auquel il est possible de télécharger les données correspond en France aux anciennes régions. Pour travailler sur Compiègne, il nous faudrait d’abord télécharger l’ensemble des données OSM de la Picardie, disponibles ici, et pesant plus de 240 Mb.
Le code ci-dessous effectue les mêmes manipulations que précédemment.
Une fois les points d’intérêt récupérés, il est possible de les agréger dans une maille.
grid <-st_make_grid(res$zone, cellsize =300, square =FALSE)grid <-st_intersection(grid, res$zone)grid <-st_sf(ID =1:length(grid), geom = grid)# compterinter <-st_intersects(grid, poi_reproj, sparse =TRUE)grid$n_poi <-lengths(inter)om_map(res, title ="Autour de la mairie")mf_map(grid, var ="n_poi", type ="choro", breaks =c(0,1,3,5,10,15),alpha = .6, border =NA, leg_title ="Nombre \nd'aménités",leg_val_rnd =0, leg_pos ="topright", add =TRUE)mf_title("Aménités gustatives - données carroyées")
6. Lissage (KDE)
Le package spatstat(Baddeley, Turner, et Rubak 2025) regroupe de nombreuses fonctions pour analyser des semis de points. Il génère des objets de classe ppp(planar point pattern), dans des fenêtres d’observation, de classe owin(observation window).
La fonction density.ppp() calcule des valeurs à partir d’une estimation par noyau (KDE). Elle prend en entrée une portée (sigma), ainsi qu’une résolution au sein de laquelle sera calculée la densité, en pixels (eps).
p <-as.ppp(X =st_coordinates(poi_reproj), W =as.owin(res$zone))ds <-density.ppp(x = p, sigma =100, eps =10, positive =TRUE)
Une fois les densités calculées, nous créons un objet de type SpatRaster. La densité étant calculée sur un mètre carré, nous la multiplions par 10.000 afin de récupérer une densité par hectare.
# Calcul densité de restaurants par hectarer <-rast(ds) *100*100crs(r) <-st_crs(poi_reproj)$wktr
class : SpatRaster
dimensions : 400, 400, 1 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
extent : 312629.1, 316629.1, 6344049, 6348049 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857)
source(s) : memory
name : lyr.1
min value : 2.225074e-304
max value : 1.970817e+00
mf_map(res$zone, col ="#f2efe9", border =NA)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5,add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(r[[2]], var ="lyr.1", type ="choro", border =NA, breaks = bks,leg_title ="Densité\n (n/ha)", pal = pal, leg_pos ="bottomleft",add =TRUE)mf_map(poi_reproj[poi_reproj$amenity =="restaurant" ,],pch =20, cex = .2, col ="black", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Densité de restaurants")
Effet de bord
7. Temps d’accès
Le package osrm(Giraud 2024) interface l’engin de routage OSRM et propose différentes fonctions pour le calcul de temps de trajet et d’itinéraires. Ces calculs s’effectuent sur la base de profils .lua qui pour chaque mode de transport définissent des algorithmes permettant la sélection de telle ou telle route, une vitesse par défaut selon le type de route, etc. Le profil voiture est disponible ici.
D’après la documentation de l’API OSRM, le nombre de requêtes sur le serveur de démo est limité à 512 par connection à l’API, espacées d’au moins 5 secondes. Pour des calculs plus volumineux (échelle nationale voire européenne), il est possible d’installer une instance OSRM sur son propre serveur (voir ici).
7.1. Depuis la mairie
Pour calculer les temps de trajet (en minutes) et les distances (en mètres) entre la mairie de Compiègne et l’ensemble des restaurants, bars et cafés situés dans un rayon de 2000m, nous utilisons la fonction osrmTable(), qui présente en sortie une liste composée de deux dataframes correspondant aux coordonnées des points d’origine et de destination, et d’un vecteur correspondant à la métrique calculée. Il est possible de calculer les distances et temps de trajet entre plusieurs couples origine / destination en une seule requête.
Nous construisons d’abord une matrice O/D composées de coordonnées X/Y des points d’origine et de destination.
# Création de la matrice O/Do <-data.frame(X =st_coordinates(pt)[, 1],Y =st_coordinates(pt)[, 2])row.names(o) <-1:nrow(o)d <-data.frame(X =st_coordinates(poi)[, 1],Y =st_coordinates(poi)[, 2])# 1 requête : 1 point d'origine, 92 points de destinationmat_dist <-osrmTable(src = o,dst = d,measure ="distance",osrm.profile ='foot')# Intervalle de 5 secondes entre les requêtesSys.sleep(5)mat_dur <-osrmTable(src = o,dst = d,measure ="duration",osrm.profile ='foot')dist <-t(mat_dist$distances)dur <-t(mat_dur$durations)poi$dist <- distpoi$dur <- durapply(poi$dur, 2, summary)
1
Min. 0.200000
1st Qu. 2.750000
Median 3.950000
Mean 6.087805
3rd Qu. 8.050000
Max. 22.700000
apply(poi$dist, 2, summary)
1
Min. 16.0000
1st Qu. 204.0000
Median 295.0000
Mean 453.6098
3rd Qu. 596.0000
Max. 1693.0000
Jouer la requête ou la passer en eval = FALSE ?
Montrer l’output de osrmTable() ?
Il est possible d’atteindre 41 restaurants, cafés ou bars en moins de 4 min à pieds depuis la mairie de Compiègne, et les trois quarts se situent à moins de 600m à pieds.
7.2. Depuis la grille
Afin de construire des indicateurs d’accessibilité, il est nécessaire de calculer les temps de trajet vers tous les restaurants / bars / cafés de la ville depuis chaque point de grille.
nrow(grid) *nrow(d)
[1] 15498
Nous avons en tout 15498 couples origine / destination. Pour éviter de surcharger le serveur de démo, nous découpons nos requêtes en autant de destinations. Ainsi, chaque requête calculera le temps de trajet entre nos 189 points d’origine et un point de destination.
Le code pour effectuer les calculs est présenté ci-dessous, mais a été joués en local. Les résultats sont enregistrés dans le dossier data.
osrmNearest()
Pour que l’engin de routage calcule un itinéraire, les points d’origine et de destination doivent être localisés sur une route. Si le point d’entrée ne se situe par sur une route, il sera déplacé sur le point le plus proche situé sur une route, en distance euclidienne. Cela peut avoir un impact sur l’itinéraire final ainsi que les temps de trajet, plus ou moins important selon le type d’espace (urbain ou rural).
# Coordonnées des points d'origine : centroïdes des cellules de la grilleo <-st_centroid(grid) |>st_transform(crs =4326)o$X <-st_coordinates(o)[, 1]o$Y <-st_coordinates(o)[, 2]o <-st_set_geometry(o[, c("X", "Y")], NULL)
Nous pouvons maintenant calculer plusieurs indicateurs d’accessibilité : temps de trajet au restaurant le plus proche, au deuxième, au troisième (question du choix), nombre de restaurants à moins de 5 minutes, 10 minutes, 15 minutes, etc.
mat <-read.csv("data/matfoot.csv")# Restaurant le plus proche pour chaque point de grillepoi_min <-aggregate(dur_foot ~ ID_ORIG, mat, FUN = min, na.rm =TRUE)# Récupérer l'identifiant de destination associépoi_min <-do.call(rbind,lapply(split(mat, mat$ID_ORIG),function(x) x[which.min(x$dur_foot), ]))# Cartographierpoi_min_sf <-merge(grid, poi_min,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)om_map(res, title ="Temps d'accès au restaurant le plus proche")mf_map(poi_min_sf, var ="dur_foot", type ="choro", breaks =seq(0,15,3),alpha = .6, border =NA, leg_title ="Temps\n(min.)",leg_val_rnd =0, pal ="Greens", leg_pos ="topright", add =TRUE)
Valhalla
Valhalla est un engin de routage open source basé sur les données OSM, au même titre que OSRM. Il fonctionne sur un système de tuiles contenant différents niveaux de routes, ainsi que les graphes routiers associés (les tuiles de niveau 0 contiennent les autoroutes, celles de niveau 2 les routes départementales). Les routes sont calculées au moment de la requête (tandis qu’elles sont pré-calculées sur OSRM) via différents algorithmes (A* pour les plus courts chemins, Contraction Hierarchies pour les longues distances). Valhalla propose de nombreux profils (vélo de location, multimodalité, taxi), ainsi que la possibilité de prendre en compte les pentes.
# Nombre de poi à moins de 5 minutesmat5 <- mat[mat$dur_foot <5 ,]mat5$n <-1mat5 <-aggregate(list(n_resto = mat5$n),list(ID_ORIG = mat5$ID_ORIG),FUN = sum)grid <-merge(grid[, "ID"], mat5,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)grid$n_resto[is.na(grid$n_resto)] <-0om_map(res, title ="Restaurants à moins de 5 minutes")mf_map(grid, var ="n_resto", type ="choro", breaks =c(0,1,5,10,30,40),alpha = .6, border =NA, leg_title ="Nb de\nrestaurants",leg_val_rnd =0, pal ="Reds", leg_pos ="topright", add =TRUE)
Le package osrm propose aussi des isochrones, à l’aide de la fonction osrmIsochrone().
Reading layer `isos' from data source
`C:\github\osm-elementr-2025-application\data\mat.gpkg' using driver `GPKG'
Simple feature collection with 6 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2.798778 ymin: 49.39991 xmax: 2.853143 ymax: 49.43269
Geodetic CRS: WGS 84
isos <-st_transform(isos, crs ="EPSG:3857")isos <-st_intersection(isos, res$zone)om_map(res, title ="Isochrones autour de la mairie de Compiègne")mf_map(isos, type ="typo", var ="isomax", alpha = .6, pal ="Geyser",border ="white", leg_pos =NA, add =TRUE)
8. Itinéraires
Pour aller plus loin, il est possible de récupérer les tracés des itinéraires, à l’aide de la fonction osrmRoute(). Attention, l’API route d’OSRM permet de calculer des itinéraires seulement entre deux points. Une requête correspond donc à un couple origine / destination.
Le package stplanr(Lovelace, Ellison, et Morgan 2024) permet initialement de modéliser et de travailler sur des systèmes de transport. Nous l’utilisons ici à des fins cartographiques, pour découper des lignes en segments.
routes <-st_read("data/mat.gpkg", layer ="routes", quiet =TRUE)routes <-st_transform(routes, crs ="EPSG:3857")routes$n <-1routes$distance <- routes$distance *1000# La fonction overline permet de compter le nombre d'occurence des tronçonsseg <-overline(routes, attrib ="n", fun = sum)routes <- routes[order(routes$distance, decreasing =TRUE), ]l_seg <-list()# La fonction line_segment découpe un objet sf LINESTRING en segments réguliersfor(i in1:nrow(routes)){ l_seg[[i]] <-line_segment( routes[i, ],segment_length =max(routes$distance) /20)}pal <-mf_get_pal(nrow(l_seg[[1]]), rev =TRUE, palette ="Blues 3")om_map(res, title ="Routes", theme ="dark")for(i innrow(routes):1){mf_map(l_seg[[i]], col = pal, lwd =1.6, add =TRUE)}leg(type ="cont", val =c(min(routes$distance), 1000, max(routes$distance)),pal = pal, pos ="left", title ="Distance (m)")
Bonus : la tournée des pizzeria
Je me trouve à la mairie de Compiègne et souhaite goûter toutes les pizzas de la ville en un minimum de temps. Avec la fonction osrmTrip(), c’est possible !
Cette fonction répond au problème du voyageur de commerce, qui consiste à déterminer le plus court circuit passant par chaque point une seule fois.
# Coordonnées des pizzeriapizz <- poi[!is.na(poi$cuisine) & poi$cuisine =="pizza" , "geometry"]# Ajouter le point de départ : la mairiepizz <-rbind(pt[1, "geometry"], pizz)pizz$id <-1:nrow(pizz)trip <-osrmTrip(loc = pizz,overview ="full",osrm.profile ="foot")trip <- trip[[1]]$tripst_write(trip, "data/mat.gpkg", layer ="trip")
trip <-st_read("data/mat.gpkg", layer ="trip", quiet =TRUE)trips <-st_transform(trip, crs ="EPSG:3857")# Pour les labels : récupérer le début de la lignestart <- lwgeom::st_startpoint(trips)start <-do.call(rbind, start)start <-data.frame(id =rownames(trips), X = start[,1], Y = start[,2])start <-st_as_sf(start, coords =c("X", "Y"), crs ="EPSG:3857")pizz <- poi[!is.na(poi$cuisine) & poi$cuisine =="pizza" ,]|>st_transform(crs ="EPSG:3857")mf_map(res$urban, col ="#a83c0a", border ="#e0dfdf", lwd = .5)mf_map(res$green, col ="#569128", border ="#569128", lwd = .5, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey80", lty =2, lwd = .2, add =TRUE)mf_map(res$building, col ="#942222", border ="#942222", lwd = .5, add=TRUE)mf_map(res$road, col ="grey80", border =NA, lwd = .2, add =TRUE)mf_map(res$street, col ="grey80", border =NA, lwd = .2, add =TRUE)mf_map(res$zone, col ="#a83c0a33", alpha = .1, border =NA, add =TRUE)mf_map(res$zone, col =NA, border ="#d7b578", lwd =15, add =TRUE) mf_map(trips, col ="black", lwd =3, lty =3, add =TRUE)mf_map(start[start$id !=1, ], pch =20, cex =2, col ="black", add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_label(start, var ="id", cex =1.2, overlap =TRUE, pos =4, col ="black", halo =TRUE)mf_title("La route de la pizza")mf_arrow()mf_scale(size =500, scale_units ="m")mf_credits("OpenStreetMap contributors, 2025", pos ="bottomleft")
sum(trip$duration)
[1] 67.045
D’après les données OSM, à peine plus d’une heure suffisent pour faire la tournée des pizzas à Compiègne.
Cambon, Jesse, Diego Hernangómez, Christopher Belanger, Daniel Possenriede, et Otto Hansen. 2025. tidygeocoder: An intuitive interface for getting data from geocoding services.https://CRAN.R-project.org/package=tidygeocoder.
Padgham, Mark, Bob Rudis, Robin Lovelace, Maëlle Salmon, et Joan Maspons. 2023. osmdata: Import ’OpenStreetMap’ Data as Simple Features or Spatial Objects. https://CRAN.R-project.org/package=osmdata.
Code source
---title: "Utiliser OpenStreetMap avec R - mise en pratique"author: - name: Louis Laurian orcid: 0009-0002-7275-5536 affiliations: - ref: riate - name: Ronan Ysebaert orcid: 0000-0002-7344-5911 - name: Timothée Giraud orcid: 0000-0002-1932-3323affiliations: - id: riate name: Centre pour l’analyse spatiale et la géovisualisation (UAR RIATE) city: Paris state: France url: https://riate.cnrs.fr/date: 13 05 2025date-format: "D MMM YYYY"lang: frformat: html: embed-resources: true smooth-scroll: true fontsize: 0.9em code-tools: true code-fold: false toc: true toc-title: Résumé toc-depth: 3 toc-location: right bibliography: bib.bib linkcolor: "#d31329" css: custom.csseditor_options: chunk_output_type: console---## ObjectifsCe document présente un usage possible des packages R qui interfacent OpenStreetMap (OSM) : **géocoder une adresse**, **import de tuiles et de couches géographiques**, **import de points d'intérêt**, **routage** et **calcul d'itinéraires**.Il propose aussi des **méthodes spatiales** (carroyages, KDE, segments de lignes) et des **visualisations associées** pour restituer graphiquement la structure spatiale de semis de points ou de lignes. En **utilisant exclusivement des informations issues d'OSM**, l'exercice consiste à **identifier les secteurs accessibles à pieds dans le voisinage géographique d'un lieu défini par une adresse**. ## Reproduire l'analysePour reproduire l'analyse, télécharger l'archive ZIP associée au dépôt [https://github.com/riatecom/osm-elementr-2025-application](https://github.com/riatecom/osm-elementr-2025-application), ouvrir le fichier proj.Rproj puis le fichier script.R qui reproduit en tout point le contenu du support de formation. Cette analyse est reproductible en n'importe quel lieu, avec une adresse (adr) et un nom de ville (lib). Néanmoins les temps de calcul utilisant le serveur de démonstration d'OSRM ont été paramétrés pour ne pas le surcharger. L'exécution de ces blocs de code prennent du temps. En conséquence, 6 jeux de données ont été préparés en amont à partir de l'adresse des gares SNCF de Lille, Toulouse, Grenoble, Nantes, Clermont-Ferrand et Compiègne, ainsi que de la station de métro Front Populaire. Les blocs de code relevant des calculs d'itinéraires / temps de trajets sont placés en commentaires dans le fichier script.R, mais sont fonctionnels.```{r}cs <-read.csv("data/case_studies.csv", encoding ="UTF-8", sep =",")cs```Nous présentons ici les résultats depuis la **gare SNCF de Compiègne**. ```{r}cs <- cs[cs$lib =="Compiegne",]adr <- cs$adrlib <- cs$lib```## 0. Packages```{r, warning = FALSE, message = FALSE}# 1. Packages qui interfacent des données OSM# install.packages('maposm', repos = 'https://riatelab.r-universe.dev')library(tidygeocoder) # géocodagelibrary(maptiles) # import de tuiles OSM (raster)library(osmdata) # import de données OSM (vecteur)library(maposm) # import de données OSM (couches géo)library(osrm) # calcul d'itinéraires# 2. Autres packages utilitaireslibrary(sf) # manipulation de données vectorielleslibrary(terra) # manipulation de données rasterlibrary(mapsf) # cartographie thématiquelibrary(maplegend) # légendeslibrary(spatstat) # analyse de semis de pointslibrary(stplanr) # segmentiser des ligneslibrary(mapview) # cartographie interactivelibrary(lwgeom)```## 1. GéocodageLe package `tidygeocoder`[@R-tidygeocoder] permet d'utiliser des services de géocodage en ligne. Par défaut, le package utilise [Nominatim](https://nominatim.org/) qui se base sur les données OSM.```{r}loc <-"Mairie de Compiègne"pt <-data.frame(address ="28 Place de l'Hôtel de Ville, 60200 Compiègne")pt <-geocode(.tbl = pt, address ="address", quiet =TRUE) |>st_as_sf(coords =c("long", "lat"), crs =4326)pt```## 2. Import de tuiles OSMLe package `maptiles`[@R-maptiles] permet de télécharger des tuiles chez différents fournisseurs, dont OpenStreetMap, au **format raster.** [20 niveaux de zoom](https://wiki.openstreetmap.org/wiki/Zoom_levels) sont disponibles, selon la taille de l'objet géographique à cartographier. Un zoom de niveau 10 correspond à une échelle 1/500.000.Plusieurs modèles de tuiles sont disponibles [ici](https://wiki.openstreetmap.org/wiki/Raster_tile_providers), certains nécessitent une inscription.```{r}emprise <-st_buffer(pt, 3000) |>st_bbox()tiles <-get_tiles(x = emprise,project =FALSE,crop =TRUE,cachedir ="cache",zoom =14)pt_reproj <-st_transform(pt, crs ="EPSG:3857")mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt = loc)mf_title("Tuiles OSM - Défaut")# Autre exemple de fournisseurtiles <-get_tiles(x = emprise,project =FALSE,crop =TRUE,provider ="OpenStreetMap.HOT",zoom =14)mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt = loc)mf_title("Tuiles OSM - OpenStreetmap.HOT")```## 3. Import de couches géographiques OSMPour la cartographie, le package `maposm` (pas sur le CRAN) se base sur les fonctions du package `osmdata`[@R-osmdata] (partie suivante) pour télécharger les données OSM, au **format vectoriel**.La fonction principale, `om_get()`, permet de récupérer certaines couches géographiques utiles (routes, parcs, bâti, etc) autour d'un point, dans un rayon donné. Le résultat prend la forme d'une liste d'objets sf.```{r}#| fig-width: 7#| fig-height: 7res <-om_get(x =c(st_coordinates(pt)[1],st_coordinates(pt)[2]),r =2000)names(res)mf_map(res$zone, col ="#f2efe9", border =NA, add =FALSE)mf_map(res$urban, col ="#e0dfdf", border ="#e0dfdf", lwd = .5, add =TRUE)mf_map(res$green, col ="#c8facc", border ="#c8facc", lwd = .5, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$building, col ="#d9d0c9", border ="#c6bab1", lwd = .5, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Autour de la mairie")mf_arrow()mf_scale(size =500, scale_units ="m")mf_credits("OpenStreetMap contributors, 2025", pos ="bottomleft")```## 4. Export de tags avec osmdataLe package `osmdata` permet d'extraire des données OSM à partir de l'API [Overpass](https://overpass-api.de/). Il nécessite dans un premier temps d'établir une *bounding box* au sein de laquelle seront extraites les données, puis de l'exécuter (fonction `add_osm_feature()`). Le résultat est une liste, qu'il faut ensuite convertir en objet sf.Les objets OSM prennent la forme de *nodes* (points), *ways* (lignes ou polygones), ou *relations* (relation entre plusieurs objets). Ces éléments sont décrits par un ensemble de *tags*, recensées sous forme de clé-valeurs (*key-values*). Ici, nous recherchons dans la clé **amenity** les points (*nodes*) correspondant aux restaurants, bars et cafés. En dehors des clés-valeurs principales, d'autres attributs peuvent être présents, comme en l'occurence le type de cuisine, ou bien les horaires.```{r, warning = FALSE}# Définition d'une bounding box (emprise Compiègne)q <- opq(bbox = emprise, osm_types = "node")# Extraction des restaurants, bars et cafésreq <- add_osm_feature(opq = q, key = 'amenity', value = c("bar", "cafe", "restaurant"))poi <- osmdata_sf(req)poi <- poi$osm_points# Sélectionner les variablespoi <- poi[, c("osm_id", "name", "amenity", "cuisine")]# Intersectionpoi <- st_intersection(poi, st_transform(res$zone, crs = 4326))# Reprojectionpoi_reproj <- st_transform(poi , crs = "EPSG:3857")```La fonction `om_map()` affiche les couches géographiques récupérées.```{r}#| fig-width: 7#| fig-height: 7# Affichageom_map(res, title ="Autour de la mairie")mf_map(poi_reproj, type ="typo", var ="amenity", cex = .7, pch =15, leg_pos ="topright", leg_val_cex = .9, leg_title_cex = .9,add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, lwd =2, col ="darkblue",add =TRUE)```Il est possible d'utiliser le package `mapview`[@R-mapview] pour afficher les points ainsi que leurs attributs, de manière interactive.```{r}mapview(res$zone, alpha.regions =0) +mapview(poi)```::: callout-note### osmextractLe package `osmextract`[@R-osmextract] permet d’extraire des données directement depuis une base de données OSM, sans avoir à passer par une API. Le plus petit niveau géographique auquel il est possible de télécharger les données correspond en France aux anciennes régions. Pour travailler sur Compiègne, il nous faudrait d'abord télécharger l'ensemble des données OSM de la Picardie, disponibles [ici](https://download.geofabrik.de/europe/france/picardie.html), et pesant plus de 240 Mb.Le code ci-dessous effectue les mêmes manipulations que précédemment.```{r, eval = FALSE}library(osmextract)osm_pt <- oe_get(place = "Picardie", layer = "points", extra_tags = "amenity", quiet = TRUE)osm_pt <- osm_pt[osm_pt$amenity %in% c("bar", "cafe", "restaurant"), ]osm_pt <- st_transform(osm_pt, crs = "EPSG:3857")osm_pt <- st_filter(osm_pt, res$zone, .predicates = st_intersects)```:::## 5. CarroyageUne fois les points d'intérêt récupérés, il est possible de les agréger dans une maille.```{r}#| fig-width: 7#| fig-height: 7grid <-st_make_grid(res$zone, cellsize =300, square =FALSE)grid <-st_intersection(grid, res$zone)grid <-st_sf(ID =1:length(grid), geom = grid)# compterinter <-st_intersects(grid, poi_reproj, sparse =TRUE)grid$n_poi <-lengths(inter)om_map(res, title ="Autour de la mairie")mf_map(grid, var ="n_poi", type ="choro", breaks =c(0,1,3,5,10,15),alpha = .6, border =NA, leg_title ="Nombre \nd'aménités",leg_val_rnd =0, leg_pos ="topright", add =TRUE)mf_title("Aménités gustatives - données carroyées")```## 6. Lissage (KDE)Le package `spatstat`[@R-spatstat] regroupe de nombreuses fonctions pour analyser des semis de points. Il génère des objets de classe **ppp** *(planar point pattern)*, dans des fenêtres d'observation, de classe **owin** *(observation window)*.La fonction `density.ppp()` calcule des valeurs à partir d'une estimation par noyau (KDE). Elle prend en entrée une portée (*sigma*), ainsi qu'une résolution au sein de laquelle sera calculée la densité, en pixels (*eps*).```{r}#| fig-cap: "@moraga2023"#| echo: falseknitr::include_graphics(rep('img/kde.png'))``````{r}p <-as.ppp(X =st_coordinates(poi_reproj), W =as.owin(res$zone))ds <-density.ppp(x = p, sigma =100, eps =10, positive =TRUE)```Une fois les densités calculées, nous créons un objet de type *SpatRaster*. La densité étant calculée sur un mètre carré, nous la multiplions par 10.000 afin de récupérer une densité par hectare.```{r, warning = FALSE}#| fig-width: 8#| fig-height: 7# Calcul densité de restaurants par hectarer <- rast(ds) * 100 * 100crs(r) <- st_crs(poi_reproj)$wktrmf_map(res$zone, col = "#f2efe9", border = NA, expandBB = c(0,.3,0,0))mf_raster(r, leg_title = "N. restaurants / ha.", leg_pos = "left", pal = c( "#75871B", "#BAB97D", "#FAF3CE", "#FDE16A", "#F9B40E", "#E88C23", "#A25933"), add = TRUE)mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(poi_reproj, pch = 20, cex = .2, col = "black", add = TRUE)mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)mf_title("La montagne alimentaire - Compiègne")```> Couleurs ? Il est possible de convertir cette couche en format vecteur. Par défaut, la fonction `as.polygons()` arrondit les valeurs des pixels à l'entier.```{r}#| fig-width: 7#| fig-height: 7# par défautr_poly <-as.polygons(r) |>st_as_sf()# tous les dixièmesr_poly <-as.polygons(r, round =TRUE, digits =1) |>st_as_sf()mf_map(res$zone, col ="#f2efe9", border =NA)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5,add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(r_poly, var ="lyr.1", breaks =quantile(r_poly$lyr.1),border =NA, lwd = .2, type ="choro", alpha =0.6,leg_title ="Densité\n (n/ha)", add =TRUE)mf_map(poi_reproj, pch =20, cex = .2, col ="black", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)```Par aménité, discrétisation commune```{r}#| fig-width: 7#| fig-height: 7# Cafés et bars ensemblepoi_reproj$amenity[poi_reproj$amenity %in%c("bar", "cafe")] <-"bar-cafe"p <-list()ds <-list()r <-list()sel <-levels(as.factor(poi_reproj$amenity))for(i in1:length(sel)){ p[[i]] <-as.ppp(X =st_coordinates(poi_reproj[poi_reproj$amenity == sel[i] ,]),W =as.owin(res$zone)) ds[[i]] <-density.ppp(x = p[[i]], sigma =100, eps =10,positive =TRUE) r[[i]] <-rast(ds[[i]]) *100*100crs(r[[i]]) <-st_crs(poi_reproj)$wkt r[[i]] <-as.polygons(r[[i]], round =TRUE, digits =1) |>st_as_sf()}# Récupérer l'amplitudeamp <-do.call(rbind, r)amp <-sort(c(unique(amp$lyr.1)))bks <-quantile(amp)pal <-mf_get_pal(n =length(bks) -1, palette ="Reds", alpha = .6, rev =TRUE)mf_map(res$zone, col ="#f2efe9", border =NA)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5,add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(r[[1]], var ="lyr.1", breaks = bks, type ="choro", border =NA,leg_title ="Densité\n (n/ha)", leg_pos ="bottomleft", pal = pal, add =TRUE)mf_map(poi_reproj[poi_reproj$amenity =="bar-cafe" ,],pch =20, cex = .2, col ="black", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Densité de bars et cafés")mf_map(res$zone, col ="#f2efe9", border =NA)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5,add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(r[[2]], var ="lyr.1", type ="choro", border =NA, breaks = bks,leg_title ="Densité\n (n/ha)", pal = pal, leg_pos ="bottomleft",add =TRUE)mf_map(poi_reproj[poi_reproj$amenity =="restaurant" ,],pch =20, cex = .2, col ="black", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Densité de restaurants")```> Effet de bord## 7. Temps d'accèsLe package `osrm`[@R-osrm] interface l'engin de routage [OSRM](https://project-osrm.org/) et propose différentes fonctions pour le calcul de temps de trajet et d'itinéraires. Ces calculs s'effectuent sur la base de **profils .lua** qui pour chaque mode de transport définissent des algorithmes permettant la sélection de telle ou telle route, une vitesse par défaut selon le type de route, etc. Le profil voiture est disponible [ici](https://github.com/Project-OSRM/osrm-backend/blob/master/profiles/car.lua).D'après la [documentation de l'API OSRM](https://project-osrm.org/docs/v5.24.0/api/#), le nombre de requêtes **sur le serveur de démo** est limité à **512** par connection à l'API, espacées d'au moins **5 secondes**.Pour des calculs plus volumineux (échelle nationale voire européenne), il est possible d'installer une instance OSRM sur son propre serveur (voir [ici](https://rcarto.github.io/posts/build_osrm_server/)).### 7.1. Depuis la mairiePour calculer les temps de trajet (en minutes) et les distances (en mètres) entre la mairie de Compiègne et l'ensemble des restaurants, bars et cafés situés dans un rayon de 2000m, nous utilisons la fonction `osrmTable()`, qui présente en sortie une liste composée de deux dataframes correspondant aux coordonnées des points d'origine et de destination, et d'un vecteur correspondant à la métrique calculée. **Il est possible de calculer les distances et temps de trajet entre plusieurs couples origine / destination en une seule requête.**Nous construisons d'abord une matrice O/D composées de coordonnées X/Y des points d'origine et de destination.```{r}# Création de la matrice O/Do <-data.frame(X =st_coordinates(pt)[, 1],Y =st_coordinates(pt)[, 2])row.names(o) <-1:nrow(o)d <-data.frame(X =st_coordinates(poi)[, 1],Y =st_coordinates(poi)[, 2])# 1 requête : 1 point d'origine, 92 points de destinationmat_dist <-osrmTable(src = o,dst = d,measure ="distance",osrm.profile ='foot')# Intervalle de 5 secondes entre les requêtesSys.sleep(5)mat_dur <-osrmTable(src = o,dst = d,measure ="duration",osrm.profile ='foot')dist <-t(mat_dist$distances)dur <-t(mat_dur$durations)poi$dist <- distpoi$dur <- durapply(poi$dur, 2, summary)apply(poi$dist, 2, summary)```> Jouer la requête ou la passer en eval = FALSE ?> Montrer l'output de osrmTable() ?Il est possible d'atteindre `r nrow(d)/2` restaurants, cafés ou bars en moins de `r round(median(poi$dur),0)` min à pieds depuis la mairie de Compiègne, et les trois quarts se situent à moins de `r round(apply(poi$dist, 2, summary)[5],-1)`m à pieds.### 7.2. Depuis la grilleAfin de construire des indicateurs d'accessibilité, il est nécessaire de calculer les temps de trajet vers tous les restaurants / bars / cafés de la ville depuis chaque point de grille. ```{r}nrow(grid) *nrow(d)```Nous avons en tout `r nrow(grid) * nrow(d)` couples origine / destination. Pour éviter de surcharger le serveur de démo, nous découpons nos requêtes en autant de destinations. Ainsi, chaque requête calculera le temps de trajet entre nos 189 points d'origine et un point de destination. Le code pour effectuer les calculs est présenté ci-dessous, mais a été joués en local. Les résultats sont enregistrés dans le dossier *data*.::: callout-important#### osrmNearest()Pour que l'engin de routage calcule un itinéraire, les points d'origine et de destination doivent être localisés sur une route. Si le point d'entrée ne se situe par sur une route, il sera déplacé sur le point le plus proche situé sur une route, en distance euclidienne. Cela peut avoir un impact sur l'itinéraire final ainsi que les temps de trajet, plus ou moins important selon le type d'espace (urbain ou rural).:::```{r, warning = FALSE}# Coordonnées des points d'origine : centroïdes des cellules de la grilleo <- st_centroid(grid) |> st_transform(crs = 4326)o$X <- st_coordinates(o)[, 1]o$Y <- st_coordinates(o)[, 2]o <- st_set_geometry(o[, c("X", "Y")], NULL)``````{r, eval = FALSE}# Les points de destination restent les mêmesmatfoot <- list()for(i in 1:nrow(d)){ matfoot[[i]] <- osrmTable(src = o[, c("X", "Y")], dst = d[i , c("X", "Y")], measure = "duration", osrm.profile = 'foot') Sys.sleep(5) matfoot[[i]]$ID_ORIG <- c(1:nrow(o)) matfoot[[i]]$ID_DEST <- rep(i, nrow(o)) matfoot[[i]] <- data.frame("ID_ORIG" = matfoot[[i]]$ID_ORIG, "ID_DEST" = matfoot[[i]]$ID_DEST, "duration" = matfoot[[i]]$durations) colnames(matfoot[[i]]) <- c("ID_ORIG", "ID_DEST", "dur_foot")}matfoot <- do.call(rbind, matfoot)write.csv(matfoot, "data/matfoot.csv", row.names = FALSE)```Nous pouvons maintenant calculer plusieurs indicateurs d'accessibilité : temps de trajet au restaurant le plus proche, au deuxième, au troisième (question du choix), nombre de restaurants à moins de 5 minutes, 10 minutes, 15 minutes, etc.```{r}#| fig-width: 7#| fig-height: 7mat <-read.csv("data/matfoot.csv")# Restaurant le plus proche pour chaque point de grillepoi_min <-aggregate(dur_foot ~ ID_ORIG, mat, FUN = min, na.rm =TRUE)# Récupérer l'identifiant de destination associépoi_min <-do.call(rbind,lapply(split(mat, mat$ID_ORIG),function(x) x[which.min(x$dur_foot), ]))# Cartographierpoi_min_sf <-merge(grid, poi_min,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)om_map(res, title ="Temps d'accès au restaurant le plus proche")mf_map(poi_min_sf, var ="dur_foot", type ="choro", breaks =seq(0,15,3),alpha = .6, border =NA, leg_title ="Temps\n(min.)",leg_val_rnd =0, pal ="Greens", leg_pos ="topright", add =TRUE)```::: callout-important#### ValhallaValhalla est un engin de routage open source basé sur les données OSM, au même titre que OSRM. Il fonctionne sur un système de **tuiles** contenant différents niveaux de routes, ainsi que les graphes routiers associés (les tuiles de niveau 0 contiennent les autoroutes, celles de niveau 2 les routes départementales). Les routes sont calculées au moment de la requête (tandis qu'elles sont pré-calculées sur OSRM) via différents algorithmes (A* pour les plus courts chemins, *Contraction Hierarchies* pour les longues distances). Valhalla propose de nombreux profils (vélo de location, multimodalité, taxi), ainsi que la possibilité de prendre en compte les pentes.:::```{r}#| fig-width: 7#| fig-height: 7# Nombre de poi à moins de 5 minutesmat5 <- mat[mat$dur_foot <5 ,]mat5$n <-1mat5 <-aggregate(list(n_resto = mat5$n),list(ID_ORIG = mat5$ID_ORIG),FUN = sum)grid <-merge(grid[, "ID"], mat5,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)grid$n_resto[is.na(grid$n_resto)] <-0om_map(res, title ="Restaurants à moins de 5 minutes")mf_map(grid, var ="n_resto", type ="choro", breaks =c(0,1,5,10,30,40),alpha = .6, border =NA, leg_title ="Nb de\nrestaurants",leg_val_rnd =0, pal ="Reds", leg_pos ="topright", add =TRUE)```Le package osrm propose aussi des isochrones, à l'aide de la fonction `osrmIsochrone()`.```{r, eval = FALSE}isos <- osrmIsochrone(loc = o, breaks = seq(0,30,5), res = 60, osrm.profile = "foot")st_write(isos, "data/mat.gpkg", layer = "isos")``````{r, warning = FALSE}#| fig-width: 7#| fig-height: 7isos <- st_read("data/mat.gpkg", layer = "isos")isos <- st_transform(isos, crs = "EPSG:3857")isos <- st_intersection(isos, res$zone)om_map(res, title = "Isochrones autour de la mairie de Compiègne")mf_map(isos, type = "typo", var = "isomax", alpha = .6, pal = "Geyser", border = "white", leg_pos = NA, add = TRUE)```## 8. ItinérairesPour aller plus loin, il est possible de récupérer les tracés des itinéraires, à l'aide de la fonction `osrmRoute()`. Attention, l'API route d'OSRM permet de calculer des itinéraires seulement entre deux points. **Une requête correspond donc à un couple origine / destination.**```{r, eval = FALSE}routes <- list()o <- data.frame(X = st_coordinates(pt)[, 1], Y = st_coordinates(pt)[, 2])row.names(o) <- 1for(i in 1:nrow(d)){ routes[[i]] <- osrmRoute(src = o, dst = d[i,], overview = "full", osrm.profile = 'foot') Sys.sleep(5)}routes <- do.call(rbind, routes)st_write(routes, "data/mat.gpkg", layer = "routes")```Le package `stplanr`[@R-stplanr] permet initialement de modéliser et de travailler sur des systèmes de transport. Nous l'utilisons ici à des fins cartographiques, pour découper des lignes en segments.```{r, warning = FALSE}#| fig-width: 8#| fig-height: 7routes <- st_read("data/mat.gpkg", layer = "routes", quiet = TRUE)routes <- st_transform(routes, crs = "EPSG:3857")routes$n <- 1routes$distance <- routes$distance * 1000# La fonction overline permet de compter le nombre d'occurence des tronçonsseg <- overline(routes, attrib = "n", fun = sum)routes <- routes[order(routes$distance, decreasing = TRUE), ]l_seg <- list()# La fonction line_segment découpe un objet sf LINESTRING en segments réguliersfor(i in 1:nrow(routes)){ l_seg[[i]] <- line_segment( routes[i, ], segment_length = max(routes$distance) / 20)}pal <- mf_get_pal(nrow(l_seg[[1]]), rev = TRUE, palette = "Blues 3")om_map(res, title = "Routes", theme = "dark")for(i in nrow(routes):1){ mf_map(l_seg[[i]], col = pal, lwd = 1.6, add = TRUE)}leg(type = "cont", val = c(min(routes$distance), 1000, max(routes$distance)), pal = pal, pos = "left", title = "Distance (m)")```## Bonus : la tournée des pizzeriaJe me trouve à la mairie de Compiègne et souhaite goûter toutes les pizzas de la ville en un minimum de temps. Avec la fonction `osrmTrip()`, c'est possible !Cette fonction répond au [problème du voyageur de commerce](https://fr.wikipedia.org/wiki/Probl%C3%A8me_du_voyageur_de_commerce), qui consiste à déterminer le plus court circuit passant par chaque point une seule fois.```{r, eval = FALSE}# Coordonnées des pizzeriapizz <- poi[!is.na(poi$cuisine) & poi$cuisine == "pizza" , "geometry"]# Ajouter le point de départ : la mairiepizz <- rbind(pt[1, "geometry"], pizz)pizz$id <- 1:nrow(pizz)trip <- osrmTrip(loc = pizz, overview = "full", osrm.profile = "foot")trip <- trip[[1]]$tripst_write(trip, "data/mat.gpkg", layer = "trip")``````{r}#| fig-width: 7#| fig-height: 7trip <-st_read("data/mat.gpkg", layer ="trip", quiet =TRUE)trips <-st_transform(trip, crs ="EPSG:3857")# Pour les labels : récupérer le début de la lignestart <- lwgeom::st_startpoint(trips)start <-do.call(rbind, start)start <-data.frame(id =rownames(trips), X = start[,1], Y = start[,2])start <-st_as_sf(start, coords =c("X", "Y"), crs ="EPSG:3857")pizz <- poi[!is.na(poi$cuisine) & poi$cuisine =="pizza" ,]|>st_transform(crs ="EPSG:3857")mf_map(res$urban, col ="#a83c0a", border ="#e0dfdf", lwd = .5)mf_map(res$green, col ="#569128", border ="#569128", lwd = .5, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey80", lty =2, lwd = .2, add =TRUE)mf_map(res$building, col ="#942222", border ="#942222", lwd = .5, add=TRUE)mf_map(res$road, col ="grey80", border =NA, lwd = .2, add =TRUE)mf_map(res$street, col ="grey80", border =NA, lwd = .2, add =TRUE)mf_map(res$zone, col ="#a83c0a33", alpha = .1, border =NA, add =TRUE)mf_map(res$zone, col =NA, border ="#d7b578", lwd =15, add =TRUE) mf_map(trips, col ="black", lwd =3, lty =3, add =TRUE)mf_map(start[start$id !=1, ], pch =20, cex =2, col ="black", add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_label(start, var ="id", cex =1.2, overlap =TRUE, pos =4, col ="black", halo =TRUE)mf_title("La route de la pizza")mf_arrow()mf_scale(size =500, scale_units ="m")mf_credits("OpenStreetMap contributors, 2025", pos ="bottomleft")sum(trip$duration)```D'après les données OSM, à peine plus d'une heure suffisent pour faire la tournée des pizzas à Compiègne.
